home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / msdos / compresn / dvpeg / src / jrevdct.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-24  |  7.4 KB  |  226 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on Appendix A.2 of the book
  11.  * "Discrete Cosine Transform---Algorithms, Advantages, Applications"
  12.  * by K.R. Rao and P. Yip  (Academic Press, Inc, London, 1990).
  13.  * It uses scaled fixed-point arithmetic instead of floating point.
  14.  */
  15.  
  16. #include "jinclude.h"
  17.  
  18. /*
  19.  * This routine is specialized to the case DCTSIZE = 8.
  20.  */
  21.  
  22. #if DCTSIZE != 8
  23.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  24. #endif
  25.  
  26.  
  27. /* The poop on this scaling stuff is as follows:
  28.  *
  29.  * We have to do addition and subtraction of the integer inputs, which
  30.  * is no problem, and multiplication by fractional constants, which is
  31.  * a problem to do in integer arithmetic.  We multiply all the constants
  32.  * by DCT_SCALE and convert them to integer constants (thus retaining
  33.  * LG2_DCT_SCALE bits of precision in the constants).  After doing a
  34.  * multiplication we have to divide the product by DCT_SCALE, with proper
  35.  * rounding, to produce the correct output.  The division can be implemented
  36.  * cheaply as a right shift of LG2_DCT_SCALE bits.  The DCT equations also
  37.  * specify an additional division by 2 on the final outputs; this can be
  38.  * folded into the right-shift by shifting one more bit (see UNFIXH).
  39.  *
  40.  * If you are planning to recode this in assembler, you might want to set
  41.  * LG2_DCT_SCALE to 15.  This loses a bit of precision, but then all the
  42.  * multiplications are between 16-bit quantities (given 8-bit JSAMPLEs!)
  43.  * so you could use a signed 16x16=>32 bit multiply instruction instead of
  44.  * full 32x32 multiply.  Unfortunately there's no way to describe such a
  45.  * multiply portably in C, so we've gone for the extra bit of accuracy here.
  46.  */
  47.  
  48. #ifdef EIGHT_BIT_SAMPLES
  49. #define LG2_DCT_SCALE 16
  50. #else
  51. #define LG2_DCT_SCALE 15    /* lose a little precision to avoid overflow */
  52. #endif
  53.  
  54. #define ONE    ((INT32) 1)
  55.  
  56. #define DCT_SCALE (ONE << LG2_DCT_SCALE)
  57.  
  58. /* In some places we shift the inputs left by a couple more bits, */
  59. /* so that they can be added to fractional results without too much */
  60. /* loss of precision. */
  61. #define LG2_OVERSCALE 2
  62. #define OVERSCALE  (ONE << LG2_OVERSCALE)
  63. #define OVERSHIFT(x)  ((x) <<= LG2_OVERSCALE)
  64.  
  65. /* Scale a fractional constant by DCT_SCALE */
  66. #define FIX(x)    ((INT32) ((x) * DCT_SCALE + 0.5))
  67.  
  68. /* Scale a fractional constant by DCT_SCALE/OVERSCALE */
  69. /* Such a constant can be multiplied with an overscaled input */
  70. /* to produce something that's scaled by DCT_SCALE */
  71. #define FIXO(x)  ((INT32) ((x) * DCT_SCALE / OVERSCALE + 0.5))
  72.  
  73. /* Descale and correctly round a value that's scaled by DCT_SCALE */
  74. #define UNFIX(x)   RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1)), LG2_DCT_SCALE)
  75.  
  76. /* Same with an additional division by 2, ie, correctly rounded UNFIX(x/2) */
  77. #define UNFIXH(x)  RIGHT_SHIFT((x) + (ONE << LG2_DCT_SCALE), LG2_DCT_SCALE+1)
  78.  
  79. /* Take a value scaled by DCT_SCALE and round to integer scaled by OVERSCALE */
  80. #define UNFIXO(x)  RIGHT_SHIFT((x) + (ONE << (LG2_DCT_SCALE-1-LG2_OVERSCALE)),\
  81.                    LG2_DCT_SCALE-LG2_OVERSCALE)
  82.  
  83. /* Here are the constants we need */
  84. /* SIN_i_j is sine of i*pi/j, scaled by DCT_SCALE */
  85. /* COS_i_j is cosine of i*pi/j, scaled by DCT_SCALE */
  86.  
  87. #define SIN_1_4 FIX(0.707106781)
  88. #define COS_1_4 SIN_1_4
  89.  
  90. #define SIN_1_8 FIX(0.382683432)
  91. #define COS_1_8 FIX(0.923879533)
  92. #define SIN_3_8 COS_1_8
  93. #define COS_3_8 SIN_1_8
  94.  
  95. #define SIN_1_16 FIX(0.195090322)
  96. #define COS_1_16 FIX(0.980785280)
  97. #define SIN_7_16 COS_1_16
  98. #define COS_7_16 SIN_1_16
  99.  
  100. #define SIN_3_16 FIX(0.555570233)
  101. #define COS_3_16 FIX(0.831469612)
  102. #define SIN_5_16 COS_3_16
  103. #define COS_5_16 SIN_3_16
  104.  
  105. /* OSIN_i_j is sine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  106. /* OCOS_i_j is cosine of i*pi/j, scaled by DCT_SCALE/OVERSCALE */
  107.  
  108. #define OSIN_1_4 FIXO(0.707106781)
  109. #define OCOS_1_4 OSIN_1_4
  110.  
  111. #define OSIN_1_8 FIXO(0.382683432)
  112. #define OCOS_1_8 FIXO(0.923879533)
  113. #define OSIN_3_8 OCOS_1_8
  114. #define OCOS_3_8 OSIN_1_8
  115.  
  116. #define OSIN_1_16 FIXO(0.195090322)
  117. #define OCOS_1_16 FIXO(0.980785280)
  118. #define OSIN_7_16 OCOS_1_16
  119. #define OCOS_7_16 OSIN_1_16
  120.  
  121. #define OSIN_3_16 FIXO(0.555570233)
  122. #define OCOS_3_16 FIXO(0.831469612)
  123. #define OSIN_5_16 OCOS_3_16
  124. #define OCOS_5_16 OSIN_3_16
  125.  
  126.  
  127. /*
  128.  * Perform the inverse DCT on one block of coefficients.
  129.  *
  130.  * A 2-D IDCT can be done by 1-D IDCT on each row
  131.  * followed by 1-D IDCT on each column.
  132.  */
  133.  
  134. GLOBAL void
  135. j_rev_dct (DCTBLOCK data)
  136. {
  137.   int pass, rowctr;
  138.   register DCTELEM *inptr, *outptr;
  139.   DCTBLOCK workspace;
  140.  
  141.   /* Each iteration of the inner loop performs one 8-point 1-D IDCT.
  142.    * It reads from a *row* of the input matrix and stores into a *column*
  143.    * of the output matrix.  In the first pass, we read from the data[] array
  144.    * and store into the local workspace[].  In the second pass, we read from
  145.    * the workspace[] array and store into data[], thus performing the
  146.    * equivalent of a columnar IDCT pass with no variable array indexing.
  147.    */
  148.  
  149.   inptr = data;            /* initialize pointers for first pass */
  150.   outptr = workspace;
  151.   for (pass = 1; pass >= 0; pass--) {
  152.     for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  153.       /* many tmps have nonoverlapping lifetime -- flashy register colourers
  154.        * should be able to do this lot very well
  155.        */
  156.       INT32 in0, in1, in2, in3, in4, in5, in6, in7;
  157.       INT32 tmp10, tmp11, tmp12, tmp13;
  158.       INT32 tmp20, tmp21, tmp22, tmp23;
  159.       INT32 tmp30, tmp31;
  160.       INT32 tmp40, tmp41, tmp42, tmp43;
  161.       INT32 tmp50, tmp51, tmp52, tmp53;
  162.       SHIFT_TEMPS
  163.     
  164.         in2 = inptr[2];
  165.         in4 = inptr[4];
  166.         in5 = inptr[5];
  167.         in6 = inptr[6];
  168.         in0 = inptr[0];
  169.  
  170.         /* These values are scaled by DCT_SCALE */
  171.  
  172.         tmp10 = (in0 + in4) * COS_1_4;
  173.         tmp13 = in6 * SIN_1_8 + in2 * COS_1_8;
  174.         tmp20 = tmp10 + tmp13;
  175.         tmp23 = tmp10 - tmp13;
  176.  
  177.         tmp11 = (in0 - in4) * COS_1_4;
  178.         tmp12 = in2 * SIN_1_8 - in6 * COS_1_8;
  179.         tmp21 = tmp11 + tmp12;
  180.         tmp22 = tmp11 - tmp12;
  181.  
  182.  
  183.         /* These values are scaled by OVERSCALE */
  184.  
  185.         in3 = inptr[3];
  186.         tmp30 = UNFIXO((in3 + in5) * COS_1_4);
  187.         tmp31 = UNFIXO((in3 - in5) * COS_1_4);
  188.  
  189.         in1 = inptr[1];
  190.         OVERSHIFT(in1);
  191.         in7 = inptr[7];
  192.         OVERSHIFT(in7);
  193.  
  194.         tmp41 = in7 + tmp31;
  195.         tmp43 = in7 - tmp31;
  196.         tmp40 = in1 + tmp30;
  197.         tmp42 = in1 - tmp30;
  198.  
  199.         /* And these are scaled by DCT_SCALE */
  200.  
  201.         tmp51 = tmp40 * OSIN_1_16 - tmp41 * OCOS_1_16;
  202.  
  203.         outptr[DCTSIZE*3] = (DCTELEM) UNFIXH(tmp23 + tmp51);
  204.         outptr[DCTSIZE*4] = (DCTELEM) UNFIXH(tmp23 - tmp51);
  205.  
  206.         tmp50 = tmp40 * OCOS_1_16 + tmp41 * OSIN_1_16;
  207.         outptr[        0] = (DCTELEM) UNFIXH(tmp20 + tmp50);
  208.         outptr[DCTSIZE*7] = (DCTELEM) UNFIXH(tmp20 - tmp50);
  209.  
  210.         tmp53 = tmp42 * OSIN_5_16 - tmp43 * OCOS_5_16;
  211.         outptr[DCTSIZE*6] = (DCTELEM) UNFIXH(tmp21 - tmp53);
  212.         outptr[DCTSIZE  ] = (DCTELEM) UNFIXH(tmp21 + tmp53);
  213.  
  214.         tmp52 = tmp42 * OCOS_5_16 + tmp43 * OSIN_5_16;
  215.         outptr[DCTSIZE*2] = (DCTELEM) UNFIXH(tmp22 + tmp52);
  216.         outptr[DCTSIZE*5] = (DCTELEM) UNFIXH(tmp22 - tmp52);
  217.  
  218.         inptr += DCTSIZE;        /* advance inptr to next row */
  219.         outptr++;            /* advance outptr to next column */
  220.      }
  221.      /* end of pass; in case it was pass 1, set up for pass 2 */
  222.      inptr = workspace;
  223.      outptr = data;
  224.   }
  225. }
  226.